# File:	    aux_functions.R
# Project:	Colonial legacy and foreign aid: Decomposing the colonial bias
# Authors:	Daina Chiba, Department of Government, University of Essex
#           Tobias Heinrich, Department of Political Science, University of South Carolina
# Contact:	dchiba@essex.ac.uk
# Date:     19 November, 2018

# Description: 
# This R script file defines auxiliary functions to be called by 
# script files No.2 -No.10. 



# Functions to perform decomposition analyses for three QoIs for Tobit.

calculate.qoi <- function(data, group.var, fit.0, fit.1){

  ## Extract beta, v-cov matrix, and variables
  
  coef.0 <- summary(fit.0 $ Sol)[[2]][,3]
  coef.1 <- summary(fit.1 $ Sol)[[2]][,3]
  data.x <- data[, names(coef.0)[-1]]
  x <- cbind(1, as.matrix(data.x))
  x.0 <- x[data[group.var]==0, ]
  x.1 <- x[data[group.var]==1, ]
  sigma.1 <- sqrt(median(fit.1 $ VCV[,1]) + median(fit.1 $ VCV[,2]) + median(fit.1 $ VCV[,3]))
  sigma.0 <- sqrt(median(fit.0 $ VCV[,1]) + median(fit.0 $ VCV[,2]) + median(fit.0 $ VCV[,3]))  
  
  ## QoI 1: Pr[ Y > 0 ] = Phi (XB / sigma)
  
  xb.x0.b0 <- x.0 %*% coef.0
  xb.x0.b1 <- x.0 %*% coef.1
  xb.x1.b0 <- x.1 %*% coef.0
  xb.x1.b1 <- x.1 %*% coef.1
  
  pry.x0.b0 <- mean( pnorm( xb.x0.b0 / sigma.0 ) )
  pry.x0.b1 <- mean( pnorm( xb.x0.b1 / sigma.1 ) )
  pry.x1.b0 <- mean( pnorm( xb.x1.b0 / sigma.0 ) )
  pry.x1.b1 <- mean( pnorm( xb.x1.b1 / sigma.1 ) )
  
  ## QoI 2: E[ Y | Y > 0 ] = XB + sigma * inv.mil, where inv.mil = dnorm(xb/sibma) / pnorm(xb/sigma)
  
  inv.mil.x0.b0 <- dnorm( xb.x0.b0 / sigma.0) / pnorm( xb.x0.b0 / sigma.0)
  inv.mil.x0.b1 <- dnorm( xb.x0.b1 / sigma.1) / pnorm( xb.x0.b1 / sigma.1)
  inv.mil.x1.b0 <- dnorm( xb.x1.b0 / sigma.0) / pnorm( xb.x1.b0 / sigma.0)
  inv.mil.x1.b1 <- dnorm( xb.x1.b1 / sigma.1) / pnorm( xb.x1.b1 / sigma.1)
  
  cy.x0.b0 <- mean( xb.x0.b0 + sigma.0 * inv.mil.x0.b0 )
  cy.x0.b1 <- mean( xb.x0.b1 + sigma.1 * inv.mil.x0.b1 )
  cy.x1.b0 <- mean( xb.x1.b0 + sigma.0 * inv.mil.x1.b0 )
  cy.x1.b1 <- mean( xb.x1.b1 + sigma.1 * inv.mil.x1.b1 )
  
  ## QoI 3: E[ Y ]
  
  ey.x0.b0 <- mean( pnorm( xb.x0.b0 / sigma.0 ) * (xb.x0.b0 + sigma.0 * inv.mil.x0.b0) )
  ey.x0.b1 <- mean( pnorm( xb.x0.b1 / sigma.1 ) * (xb.x0.b1 + sigma.1 * inv.mil.x0.b1) )
  ey.x1.b0 <- mean( pnorm( xb.x1.b0 / sigma.0 ) * (xb.x1.b0 + sigma.0 * inv.mil.x1.b0) )
  ey.x1.b1 <- mean( pnorm( xb.x1.b1 / sigma.1 ) * (xb.x1.b1 + sigma.1 * inv.mil.x1.b1) )
  
  out <- data.frame( pry.x0.b0 = pry.x0.b0, pry.x0.b1 = pry.x0.b1, 
                     pry.x1.b0 = pry.x1.b0, pry.x1.b1 = pry.x1.b1, 
                     cy.x0.b0 = cy.x0.b0, cy.x0.b1 = cy.x0.b1, 
                     cy.x1.b0 = cy.x1.b0, cy.x1.b1 = cy.x1.b1,
                     ey.x0.b0 = ey.x0.b0, ey.x0.b1 = ey.x0.b1, 
                     ey.x1.b0 = ey.x1.b0, ey.x1.b1 = ey.x1.b1)
  return(out)
}



# Same function but it does the calculation for those coefficients simulated 
# from MVN(MLE_beta, var-cov matrix) to approximate the posterior.

calculate.qoi.full <- function(data, group.var, fit.0, fit.1, beta.0, beta.1){

  ## Extract beta, v-cov matrix, and variables
  
  coef.0 <- summary(fit.0 $ Sol)[[2]][,3]
  coef.1 <- summary(fit.1 $ Sol)[[2]][,3]
  data.x <- data[, names(coef.0)[-1]]
  x <- cbind(1, as.matrix(data.x))
  x.0 <- x[data[group.var]==0, ]
  x.1 <- x[data[group.var]==1, ]
  sigma.1 <- sqrt(median(fit.1 $ VCV[,1]) + median(fit.1 $ VCV[,2]) + median(fit.1 $ VCV[,3]))
  sigma.0 <- sqrt(median(fit.0 $ VCV[,1]) + median(fit.0 $ VCV[,2]) + median(fit.0 $ VCV[,3]))  
  
  ## QoI 1: Pr[ Y > 0 ] = Phi (XB / sigma)
  
  xb.x0.b0 <- x.0 %*% beta.0
  xb.x0.b1 <- x.0 %*% beta.1
  xb.x1.b0 <- x.1 %*% beta.0
  xb.x1.b1 <- x.1 %*% beta.1
  
  pry.x0.b0 <- mean( pnorm( xb.x0.b0 / sigma.0 ) )
  pry.x0.b1 <- mean( pnorm( xb.x0.b1 / sigma.1 ) )
  pry.x1.b0 <- mean( pnorm( xb.x1.b0 / sigma.0 ) )
  pry.x1.b1 <- mean( pnorm( xb.x1.b1 / sigma.1 ) )
  
  ## QoI 2: E[ Y | Y > 0 ] = XB + sigma * inv.mil, where inv.mil = dnorm(xb/sibma) / pnorm(xb/sigma)
  
  inv.mil.x0.b0 <- dnorm( xb.x0.b0 / sigma.0) / pnorm( xb.x0.b0 / sigma.0)
  inv.mil.x0.b1 <- dnorm( xb.x0.b1 / sigma.1) / pnorm( xb.x0.b1 / sigma.1)
  inv.mil.x1.b0 <- dnorm( xb.x1.b0 / sigma.0) / pnorm( xb.x1.b0 / sigma.0)
  inv.mil.x1.b1 <- dnorm( xb.x1.b1 / sigma.1) / pnorm( xb.x1.b1 / sigma.1)
  
  cy.x0.b0 <- mean( xb.x0.b0 + sigma.0 * inv.mil.x0.b0 )
  cy.x0.b1 <- mean( xb.x0.b1 + sigma.1 * inv.mil.x0.b1 )
  cy.x1.b0 <- mean( xb.x1.b0 + sigma.0 * inv.mil.x1.b0 )
  cy.x1.b1 <- mean( xb.x1.b1 + sigma.1 * inv.mil.x1.b1 )
  
  ## QoI 3: E[ Y ]
  
  ey.x0.b0 <- mean( pnorm( xb.x0.b0 / sigma.0 ) * (xb.x0.b0 + sigma.0 * inv.mil.x0.b0) )
  ey.x0.b1 <- mean( pnorm( xb.x0.b1 / sigma.1 ) * (xb.x0.b1 + sigma.1 * inv.mil.x0.b1) )
  ey.x1.b0 <- mean( pnorm( xb.x1.b0 / sigma.0 ) * (xb.x1.b0 + sigma.0 * inv.mil.x1.b0) )
  ey.x1.b1 <- mean( pnorm( xb.x1.b1 / sigma.1 ) * (xb.x1.b1 + sigma.1 * inv.mil.x1.b1) )
  
  out <- data.frame( pry.x0.b0 = pry.x0.b0, pry.x0.b1 = pry.x0.b1, 
                     pry.x1.b0 = pry.x1.b0, pry.x1.b1 = pry.x1.b1, 
                     cy.x0.b0 = cy.x0.b0, cy.x0.b1 = cy.x0.b1, 
                     cy.x1.b0 = cy.x1.b0, cy.x1.b1 = cy.x1.b1,
                     ey.x0.b0 = ey.x0.b0, ey.x0.b1 = ey.x0.b1, 
                     ey.x1.b0 = ey.x1.b0, ey.x1.b1 = ey.x1.b1)
  return(out)
}

